Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 2 December 2009 (MN IM^K style file v2.2) 



Testing spectral models for stellar populations with star 
clusters: I. Methodology 

O Roberto Cid Fernandes 1 * and Rosa M. Gonzalez Delgado 2 f 

, 1 Departamento de Fisica-CFM, Universidade Federal de Santa Catarina, P.O. Box 476, 88040-900, Florianopolis, SC, Brazil 
s ! ■ 2 Instituto de Astrofisica de Andalucia (CSIC), P.O. Box 3004, 18080 Granada, Spain 

o : 

Q ' 2009 July 

<N . 

ABSTRACT 

/""N \ High resolution spectral models for simple stellar populations (SSP) developed in the 

past few years have become a standard ingredient in studies of stellar population of 
galaxies. As more such models become available, it becomes increasingly important 
C"| \ to test them. In this and a companion paper, we test a suite of publicly available 

Qlc evolutionary synthesis models using integrated optical spectra in the blue-near-UV 

range of 27 well studied star clusters from the work of Leonardi & Rose (2003) spanning 
a wide range of ages and metallicities. Most (23) of the clusters are from the Magellanic 
clouds. This paper concentrates on methodological aspects of spectral fitting. The 
| data are fitted with SSP spectral models from Vazdekis and collaborators, based on 

the MILES library. Best-fit and Bayesian estimates of age, metallicity and extinction 
are presented, and degeneracies between these parameters are mapped. We find that 
these models can match the observed spectra very well in most cases, with small 
formal uncertainties in t. Z and Ay . In some cases, the spectral fits indicate that the 
models lack a blue old population, probably associated with the horizontal branch. 
This methodology, which is mostly based on the publicly available code starlight, 
f^S ' is extended to other sets of models in Paper II, where a comparison with properties 

derived from spatially resolved data (color-magnitude diagrams) is presented. The 
global aim of these two papers is to provide guidance to users of evolutionary synthesis 
models and empirical feedback to model makers. 
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1 INTRODUCTION lution evolutionary synthesis models, which have appeared 

. . , , , , r , , ■ in the literature in the past few years (see e.g. Bruzual 2007; 

btar clusters (SCs) are ideal test beds lor the evolution „ „ , , „„„„ , r , i • \ mi i 

. „ Gonzalez Delgado 2009 and references therein), these mod- 

of stars and stellar populations. Their key properties are , ,, ,, „,,. r , , , , , , 

f , 4 \ , ers allow the fitting of observed spectra on a A-by-A ba- 

their age It) and metallicity (Z), but extinction (Ay) also ■ ■ ,■ i, -i , i ■ r ^ ■ i 

, . , , i t sls > thus incorporating all available information. Curiously, 

plays a role m determining observed properties. The classical , . ,. , , , , , ,. 

though numerous studies have used these models to disen- 
method to estimate these properties is through the compar- , . , „ , nl , , . . , . , 

f tangle the mixture of stellar populations m galaxies (e.g., 

ison of observed color-magnitude diagrams (CMDs to the „. , „ . , , £ ^ . s , , 

v ' Cid t'ernandes et al. 2007 and references therein), refatively 

predictions of stellar evolution models. For distant SCs, how- , , , ,, , . , „ , 

little has been done on the much simpler problem of appfy- 

ever, only integrated light measurements are available, and . , . , . r r r,^ mi 

, ,. . . , ing spectral synthesis to inter the properties ol bCs. the 

the traditional approach m this case is to use broad band 

colors or spectral indices to estimate their properties. Both 



most complete study in this vein has just been published 

by Koleva et al. (2008, K08), who analyzed optical spectra 
these methods have a long history m the literature (see, e.g., r^i iui i c ii a \.< ii /nnntr\ 

J _ , ' ° of Gafactic globular clusters from the bchiavon et al. (2005) 



the conference books Lamers et al. 2004; Perez et al. 2009; 
Geisler et al. 2001) 



atlas. They find that spectral fitting provides t and Z esti- 
mates in very good agreement with those obtained through a 



A third and in principle more thorough alternative has ™ m , . ,. , r , 

, , ... CMD analysis, thus validating spectral synthesis as a usefuf 

become possible with the availability of high spectral reso- , c „„ , . ., , , , r . Trt 

J ° r tool for bC work, bimilar work, but focusing on the near-lK 

range, has been recently published by Langon et al. (2009). 

* E-mail:cid@astro. ufsc.br This paper is dedicated to methodological aspects of 

f E-mail: rosa@iaa.es spectral fitting of SCs. Detailed fits of optical spectra are 
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used to estimate t, Z and Av for 27 SCs from the sam- 
ple of Leonardi & Rose (2003; hereafter LR03). A thorough 
mapping of the parameter space is performed to quantify 
the uncertainties and degeneracies in the parameters. All 
results presented here are based on fits obtained with the 
Vazdekis-MILES single stellar populations (SSP) models. 
In a companion paper (Gonzalez Delgado & Cid Fernan- 
des 2009, Paper II) the same methodology is extended to 
a suite of other publicly available high resolution evolution- 
ary synthesis models, allowing us to evaluate their pros and 
cons in a comparative fashion, as well as to quantify the un- 
certainties in SC properties resulting from variance in the 
ingredients for the analysis. Paper II also compares the SC 
properties derived from spectral fitting to those obtained 
from CMD work on the same clusters. 

The central goal of these two papers is to provide useful 
guidance to users of evolutionary synthesis models and spec- 
tral synthesis tools, as well as empirical feedback to model 
makers. Given these goals, we deliberately restrict our anal- 
ysis to models and analysis tools as found in public distribu- 
tions, with as little as possible manipulation. We also point 
out that although this study focus entirely on SCs, its un- 
derlying motivation is to evaluate the reliability of methods 
and ingredients widely disseminated in modern analysis of 
the integrated light of galaxies (see Cid Fernandes 2007 and 
references therein). SCs are the undoubtedly the best test- 
beds for this. If SC ages can be accurately recovered with 
current models and spectral synthesis methods, this lends 
confidence to applications which seek to estimate the star 
formation histories of galaxies. Similarly, if metallicities are 
also well recovered, then chemical evolution studies are war- 
ranted. In short, establishing the accuracy with which age 
and metallicity of SCs can be derived through spectral syn- 
thesis is important to understand the systematics and limi- 
tations of studies of galaxy evolution based on the analysis 
of the fossil record encoded in their spectra. 

The SC data and SSP models are described in Sections 
[2] and [3] respectively. Section [4] describes the spectral fitting 
method. Appendices IA1 and [Bl describe technical aspects of 
the fits and simulations to evaluate the effects of noise. Re- 
sults of these fits are presented in Section [5] which also il- 
lustrates and discusses the degeneracies involved. Section [5] 
applies a Bayesian formalism to provide formal estimates of 
the uncertainties in the derived parameters. Finally, Section 
[7] summarizes our main results. 



2 DATA 

Spectra for 23 SCs in the Magellanic Clouds (3 in the SMC 
and 20 in LMC) and 4 Galactic clusters (GCs) were col- 
lected by LR03 with the 1.5m telescope at the Cerro Tololo 
Inter-American Observatory, and kindly made available to 
us. The spectra cover the wavelength range 3500-4700 A 
with a final resolution FWHMinst = 5.7 A (including the 
gaussian smoothing applied by LR03). Our analysis restricts 
the spectral range to 3650-4600 A, an interval covered by all 
available SSP models, and which skips calibration problems 
towards the edges of some of the spectra. The signal-to-noise 
ratio of these spectra, measured from the rms deviation of 
Fx in the 4010-4060 A window, range from 19 to 112, with 
a sample average S/N of 54. 



In addition to the high quality of the spectra, the moti- 
vation to re-analyze these data with a spectral fitting ap- 
proach stems from the fact that their t's and Z's were 
determined through independent means, including CMDs, 
spectral indices, colors and spectroscopy of individual stars 
(LR03). These independent estimates allow a comparative 
analysis similar to that performed by K08. A further mo- 
tivation is that, compared to the systems studied by K08, 
these SCs are generally younger (from less than 0.1 to 3 
Gyr according to LR03), which makes Z estimates more 
challenging (Wolf et al. 2007). 



3 EVOLUTIONARY SYNTHESIS MODELS 

Evolutionary tracks and stellar libraries are the two most 
important ingredients in evolutionary synthesis models. As 
a result of recent work in both fronts, there are nowadays 
several sets of models available for spectral synthesis work. 
One of our major goals is to subject these models to an em- 
pirical test. This is done in Paper II, where results obtained 
with different publicly available models are compared. For 
clarity, this paper uses only one version of these models. 

We chose to work with the set of SSP mod- 
els built by Vazdekis et al. (2009), available at 
http://www.ucm.es/info/Astrof/miles/models/models.html 
The main novelty in these models is that they incorporate 
the MILES library of Sanchez-Blazquez et al. (2006). 
MILES contains about one thousand stars spanning a large 
range of stellar parameters. The spectra cover the 3540 
to 7410 A range with a resolution of 2.3 A (FWHM). 
The spectra are well calibrated in flux and wavelength, 
and represent a significant improvement with respect to 
previous libraries in the coverage of the metallicity, number 
of giant stars and other aspects. (The ELODIE library of 
Le Borgne et al. 2004 is similar in size and coverage of 
stellar parameters, but, starting at 4000 A, does not cover 
the range covered by our data.) 

This library was incorporated into the evolutionary syn- 
thesis code of Vazdekis (1999), using "Padova 2000" evolu- 
tionary tracks (Girardi et al. 2000, 2002) and a Salpeter IMF 
(Salpeter 1955). SSP spectra were computed for Nz = 6 
metallicities: Z = 0.0004, 0.001, 0.004, 0.008, 0.019 (solar) 
and 0.03, and N t = 46 ages between 0.10 and 17.78 Gyr. 
The lack of younger ages is due to the small number of hot 
(over 15000 K) stars in MILES. Few SCs in our sample are 
young enough to be affected by this limitation. 

As with other sets of models (see Paper II) the sam- 
pling in age is practically continuous, but predictions are 
only available for a few metallicities. Model-users are thus 
limited to a coarse sampling in Z. Rigorously speaking, the 
fact that evolutionary stages occur at different times for dif- 
ferent Z's makes interpolation to a finer grid physically in- 
valid, even though reasonable approximations can be ob- 
tained. Given our goal to remain as faithful as possible to 
SSP-models in their original form, available to users in gen- 
eral, we concentrate on results obtained with the original 
Nz = 6 grid. Experiments with Z-interpolated SSP spectra 
are reported in Appendix lAl 

Throughout the rest of this paper ages will be given in 
yr, and metallicities will be quoted in a log-solar scale. In 
this notation, the Z values above correspond to (rounding 



Testing spectral models for stellar populations 3 



up to the first decimal) log Z/Z Q = -1.7, -1.3, -0.7, -0.4, 
0, and +0.2. For consistency with the notation in Paper II, 
this set of models will be referred to as "VOOs" (an acronym 
for "Vazdekis + Padova 2000 + Salpeter"). 



4 SPECTRAL FITTING: METHOD 

The goal of spectral fitting is to find a model spectrum (Ma) 
which best matches the observed one (Oa), taking into ac- 
count all A's. The fits presented below were carried out with 
version 05 of the publicly available starlight^] code, bet- 
ter known as a tool to retrieve the star formation history 
of galaxies by fitting a spectrum with a mixture of iV* SSPs 
from a base spanning wide ranges of t and Z (see Cid Fer- 
nandes et al 2005, 2009; Asari et al 2007 and the starlight 
user manual for examples and references). 

A so far unused feature of this code is that it also fits 
Oa with each of the base spectra individually, thus making 
it useful to study single-population systems like SCs. This 
section describes the parameters and technical aspects of 
these fits. The level of details is justified by the fact that 
this is the first time starlight is used to fit SCs. 

A series of technical details lurk behind the deceiving 
simplicity of fitting an observed SC spectrum. Should kine- 
matical parameters be fitted when they are not relevant? 
How many degrees of freedom are involved in a fit? What 
to do when no measure of the error in Ox is available? How 
do noise affects the fits? Should one work with a discrete 
base of SSPs or a ~ continuous one obtained through in- 
terpolation? Some of these issues are discussed below, while 
others are addressed in Appendices [A] and [B] Because of 
the general goals outlined in Section [T] our default strat- 
egy to tackle these issues is to mimic inasmuch as possible 
the way starlight is used to analyze galaxy spectra, while 
at the same time avoiding manipulations of the original SSP 
models. Variations over this strategy are discussed whenever 
relevant. 

4.1 The model spectrum and its parameters 

The equation for the model SSP spectrum fitted by 
starlight is 

Ma = xji SP (t, z)lQ-° AAviq ^- q ^ (1) 

where £ is a scaling factor, q\ = Ax/Ay is the reddening 
curve, and 

ll SP it,Z) = ^p^±®G{v^) (2) 

gives the spectrum of an SSP of age t and metallicity Z nor- 
malized at Ao = 4020 A, convolved with a gaussian filter 
centered at centered at velocity and with dispersion a*. 
The Lf SP (t, Z) spectra are taken directly from the evolu- 
tionary synthesis models. 

To find the SSP which best matches a given observed 
spectrum one may construct a base of N* = Nt X Nz = 

1 www.starlight.ufsc.br 



46 x 6 = 276 SSP model spectra. For each base component 
starlight then finds the corresponding values of Ay and x 
which best match Ox- For practical reasons, the actual fits 
were carried out using 6 single Z bases, each containing all 
Nt available ages. In the end one has a table of \ 2 , Ay and 
x for each of the Nz x Nt combinations of t and Z. 

The free parameters in these fits are t, Z, Ay and x. Of 
these, t and Z are the clearly the more fundamental ones. 
Ay also has some astrophysical value, but the scaling factor 
x is just a technical parameter, with no physical relevance 
f see 14.1.2]) . Note that t and Z can only assume the values in 
the base, while Ay and x are continuous. 



4-1.1 Extinction 

Unlike the method of K08, which concentrates on fitting ab- 
sorption lines by modeling the continuum with a high order 
polynomiun, or the "Continuum-normalized" fits of Wolf et 
al. (2007), our spectral fits do take the continuum shape into 
account, so dust effects must be considered. Considering the 
difficulty in defining the continuum around the Balmer jump 
and the 4000 A break, fits to the full non-rectified spec- 
tra seem more advisable. More importantly, unless there are 
calibration problems with the data or the models, it is ob- 
viously advantageous to take the full spectral information 
into account. 

Extinction is dealt with assuming a simple foreground 
screen model (equation[T} , unrealistic for galaxies but Ok for 
SCs (at least for the ages of our clusters). SMC, LMC and 
Galactic clusters were fitted with reddening curves appropri- 
ate to each of these galaxies (Gordon et al. 2003; Cardelli, 
Clayton & Mathis 1989). 

The V-band extinction is a free parameter in our fits, 
but a priori limits can be set which may, at least in principle, 
help constraining the fits. The extinction in the LMC and 
SMC is low. Bica & Alloin (1986) give a global extinction of 
E(B - V) ^ 0.06 and 0.03 for the SCs in LMC and SMC, 
respectively, and between 0.0 and 0.09 for the GCs in our 
sample. We have also compiled LMC extinction values using 
the web tool described by Zaritsky et al. (2004) , which gives 
Ay values along the line-of-sight to stars within a search 
radius of the target coordinates. At the positions of the LMC 
SCs, we obtain Ay ranging from 0.3 to 0.7, with a large 
dispersion: a (Ay) = 0.3-0.6. These values are much larger 
than those given by Bica & Alloin (1986), but Zaritsky et al. 
explain that, given the highly non-gaussian distribution of 
Ay values, this method can only provide a rough estimation 
of the extinction. As a whole, Ay ^ 1 is a safe upper limit 
for all SCs studied here. 

We also impose the physical limit Ay ^ 0. While natu- 
ral, this choice deserves a comment. Imposing Ay ^ pre- 
vents starlight from compensating for possible flux cal- 
ibration problems in the stellar libraries. Such problems 
are known to exist with STELIB (Le Borgne et al. 2003). 
Comparing stars in common between MILES and STELIB, 
Sanchez-Blazquez et al. (2006) find the latter to be slightly 
too red. This is the reason why ~ dustless galaxies are some- 
what better fitted with slightly negative values of Ay (of 
the order of —0.1 mag) when the STELIB-based Bruzual 
& Chariot (2003) models are used, while with MILES this 
problem goes away (Gomes 2009; Cid Fernandes et al. 2009). 
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4-1.2 Scaling factor 

Since both model and observed spectra are normalized at 
Ao, the scaling factor x in equation |T} would seem an un- 
necessary parameter. This is not strictly true. First, models 
and data are normalized in slightly different ways: While 
the SSP models are normalized exactly at Ao = 4020 A, the 
observed spectrum Ox is actually normalized by the median 
flux in the 4010-4060 A window — starlight uses this trick 
to circumvent eventual problems in the Ao pixel. Second, 
after convolution with the kinematical filter G, the 7f SP 
spectra are not exactly = 1 at Ao (see equation [2}. Regard- 
less of these technicalities, forcing a M\ — O\ normal- 
ization would imply treating Ao differently from other A's, 
ultimately making the whole fit formally dependent on the 
arbitrary choice for Ao- For these reasons, the scaling factor 
must be considered a free parameter, even though one does 
not expect large departures from x = 1. 



4-.1.3 Kinematical parameters 

Formally, and cr* could also be considered free param- 
eters. However, these are not varied during the single- 
population fits. Instead, they are fixed at the values de- 
termined during the initial multi-population fit, where the 
spectrum is computed from 



M; 



ESSP 
Xj7x 



(t jt Zj) x i(r a4A v-(9A-<JA ) (3) 



3=1 



instead of equation fTJ. As expected, we always obtain v* 
within a few kms -1 of zero, since the spectra are in the 
rest-frame. Hence, should not be considered as a truly 
free parameter. The same applies to cr*, which, for the best 
fit models, is always close to the value corresponding to the 
difference between spectral resolution of the models and the 
data — the spectral resolution of the data analyzed here is 
o~inst ~ 170 kms -1 , much larger than the actual velocity 
dispersion of SCs, so that no useful kinematical informa- 
tion can be derived. We therefore do not count er* as a free 
parameter either. 

starlight allows the user to fix er* at its expected value 
(~ o-i, ls t, in our case), and this would be a valid strategy for 
the SC data analyzed here. Given our goal of mimicking 
inasmuch as possible the way that spectral synthesis is car- 
ried out with galaxy data, we opted not to fix a*. For SC 
spectra, this can be seen as conservative choice, since a* can 
be used to compensate for a metallicity mismatch (Koleva 
et al. 2009). Increasing cr* has the effect of decreasing line 
depths, which for metal lines can mimic the effect of decreas- 
ing metallicity, and vice- versa, so fitting a low Z system with 
a higher Z model and free cr* will yield a x 2 n °t as bad as 
with a fixed a*. This effect is in fact detected in our fits. We 
find that the returned 0+ values tend to increase with Z. 
Typical values span the 150 to 200 kms -1 range, bracketing 

0~inst - 

Overall, however, the impact of our choice to let cr* free 
is insignificant. Comparing SC parameters derived with fits 
with fixed and free cr* leads to typical differences of 0.07 or 
less for log t, log Z and Ay. We are clearly not yet in a posi- 
tion to aim this level of accuracy. For instance, as shown in 



Paper II, such differences are much smaller than those result- 
ing from the use of different evolutionary synthesis models. 

4.2 Error spectrum 

The figure of merit used by starlight fits is a standard 



X 



(4) 



where ex is the error in Ox. As is often the case, we do not 
have formal values for the errors. The global amplitude of 
ex is not important for the minimization of x 2 , as it does 
not change the shape of the likelihood function, and, as ex- 
plained below, \ 2 values are rescaled anyway. The shape of 
the error spectrum, which defines whether certain pixels are 
given more weight than others, is more relevant. 

Two recipes for ex were explored. In the first one we set 
ex = O.IOa, ie, a flat error spectrum with amplitude equal 
to 10% of the average flux, such that all pixels are given 
the same weight. The second recipe is ex = O.IOa, which 
gives a larger weight to absorption lines. We have verified 
that the results obtained with these two recipes are nearly 
always identical, so only results for the second recipe are 
presented. When the two recipes yield different results the 
best fit model is clearly inadequate anyway. This happens, 
for instance, when trying to fit very metal poor SCs with 
models of higher metallicity. 



4.3 x 2 re-scaling and the effective number of 
degrees of freedom 

The lack of proper errors implies that our x 2 values do not 
have a formal statistical meaning. While this does not affect 
the identification of a best model, the absolute scale of x 2 
is relevant to assess uncertainties in the derived parameters. 
In cases like this, it is common to re-scale the errors such 
that the best model has a x 2 equal to the number of degrees 
of freedom (e.g., Barth et al. 2002; Garcia-Rissmann et al. 
2005). Using a subscript "ST" to denote the x 2 value re- 
turned by starlight and "min" to indicate the best model, 
the corrected x 2 is then given by 



2 

X =Ndof- 



XST,, 



For convenience, we define 



X -x mm XST ~ XST .777 7-7) 



X 2 ■ 

A n) in 



XsT,r. 



such that changes in x 2 scale with 8 2 : 



A 2 2 

Ax =X 



= N dof 5 x 



(5) 



(6) 



(7) 



In these equations Ndof = Nx — N par , i.e., the number 
of data points minus the number of free parameters. N par 
is well defined: There are 4 free parameters {t, Z, Ay and 
a;). But what should be used for JVa? In other words: How 
many observables are we actually fitting? 

Our fits are performed with a AA = 1 A sampling, 
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so there are 951 data points between 3650 and 4600 A. 
However, the observed spectra are highly oversampled, so 
the number of actually independent data points is <C 951. 
Counting only points separated by 1 FWHMi„,i = 5.7 A 
would yield N\ — 167, but this is still an overestimate, given 
the heavy overlap in the instrumental profiles for a spectral 
distance of just 1 FWHM inst . To be safe, only points sep- 
arated by 6er; ns t = 14.5 A will be counted as independent, 
such that the overlap in instrumental profiles occurs beyond 
±3cri„ st , being therefore insignificant. This somewhat sub- 
jective but clearly conservative choice yields Nx = 65, and 
hence Ndof — 61. (These values are slight smaller in a few 
cases due to masked windows.) This recipe only affects the 
numerical interpretation of the statistical confidence associ- 
ated to a given value of A^ 2 , which can be straightforwardly 
recomputed for any other choice of Ndof- 

In the presentation of the spectral fits two other figures 
of merit are used: the rms of the R\ = Ox — M\ residual 
and the mean value of \R\\/0\, denoted by A. These nu- 
merically similar indices give an easy-to-interpret measure 
of the quality of the fit. 



5 SPECTRAL FITTING: RESULTS 

The method outlined above was applied to the 27 SC spectra 
described in Section [2] after basic pre-processing steps, such 
as resampling to AA = 1 A, and masking bad pixels. This 
section presents the results obtained. 

First we illustrate our methodology with fixed Z fits 
fSection l5.1|) and discuss the meaning of multi-SSPs results 
for SCs (Section 15. 2) . Then, in Section [5.31 we inspect how 
much the fits constrain t, Z and Ay (a task which will 
be readdressed more rigorously in Section [6j , and investi- 
gate differences in spectral fits obtained with different Z's. 
Third, we map covariances ("degeneracies") between param- 
eters (Section . For didactic purposes, the LMC clusters 
NGC 2010 and NGC 2210 are taken as examples. Section 
15.51 presents results for all SCs. 

5.1 Fixed Z analysis 

Fig. Q] shows the fit to NGC 2010 obtained with Z fixed 
at Zq. Insignificantly better (S x 2 — 0.01) fits are obtained 
with the 1.5Zq models. The left panel shows the observed, 
model, and residual spectra. The best match is obtained 
for log* = 8.20, A v = 0.00 and x = 1.03. The middle right 
panel shows the runs of Av and function of t, while the 

bottom one shows the two indicators of fit quality discussed 
at the end of Section [4.3l The largest differences between Ox 
and M\ are found in the continuum between higher order 
Balmer lines and at the bottom of the Call K line, but 
these are relatively small deviations. Overall, the spectral 
fit is very good both in the continuum and absorption lines. 

Fig. [2] shows the results for NGC 2210 and log Z/Z GJ = 
— 1.7. This is the Z which produces the smallest \ 2 among 
all 6 Z's in the VOOs models. Since \ogZ/Z Q = -1.7 is 
also the smallest one in the grid, it is possible that models 
with even lower Z would provide better fits. In this case 
differences between Ox and Mx concentrate around 4200 A. 
The best fit is achieved for log* = 9.80 and A v = 0.10. Fits 
for younger ages require a larger extinction to compensate 



for the bluer predicted colors, but even then the quality of 
the fit deteriorates quickly as one moves away from the best 
t. Notice also that Ay saturates at the imposed limits as 
one moves away from the best model. 

5.2 The meaning of multi component fits 

Accepting that SCs are single population systems, multi- 
SSP fits are in principle not of direct interest. Still, even if 
only out of curiosity, lets have a look at them. 

The top right panel in Fig. [T] shows how the light at 
Ao is spread among SSPs of different ages in the multi-SSP 
starlight fit of NGC 2010. In this particular example, the 
multi-SSP fit is identical to the single-SSP one. In other 
words, given total freedom to mix 46 different base elements, 
starlight preferred to use only one. In the Z = 1.5-Zq fits 
(not shown), the multi-SSP fit splits into two similar parts, 
corresponding to populations within 0.1 dex of the best-fit 
single SSP age. The fit is thus well focused, such that multi 
and single SSP fits are almost identical. 

The situation is very different in NGC 2210 (Fig. [2j). 
The multi-SSP fit in this case shows two dominant compo- 
nents with ages over 2 dex apart: One with logt = 10.10, 
responsible for about 80% of the light, and other at the 
youngest age logt = 8.00, accounting for the rest. Unlike 
in Fig. [T] the difference between the multi and single-SSP 
spectra (shown in the bottom of the left panel) is now no- 
ticeable. Accordingly, the multi-SSP fit yields a \ 2 which is 
27% better than the single-SSP one (a difference which is of 
only 5% for NGC 2010). 

A plausible interpretation of this "old plus a bit of 
young" population mixture, which happens in at least two 
other clusters in our sample, is that the SSP models lack 
old blue stars. This same problem was identified by K08 in 
a different sample. They found that clusters with a blue hor- 
izontal branch (HB) are better fitted adding a set of stars 
in the T off = 6000 to 20000 K range to a pure SSP. They 
also find that this ad hoc, but physically reasonable, recipe 
produces ages in better agreement with those derived from 
CMD studies. We will return to this issue in Section [7] For 
the moment, we note that the example of NGC 2210 fully 
corroborates this interpretation: The fiducial age adopted 
by LR03 is logt = 10.1, identical to that of the oldest popu- 
lation in our multi-SSP fits, but 0.3 dex older than the one 
obtained with a single-SSP fit (logt = 9.80). 

These results show that multi-SSP fits provide an use- 
ful empirical measure of the adequacy of single-SSP fits. 
Six quantities are used to report results of the multi-SSP 
fits: The mean (logt m ) and standard deviation (<r m ) of logt 
(computed directly from the population vector x using equa- 
tions 2 and 4 of Cid Fernandes et al. 2005), the age of the 
dominant population (logT m ) and its corresponding per- 
centage light fraction (/ m ), plus the mean relative residual 
(A m ), and the difference in \ 2 between the multi and single 
SSP fits (5m), defined as in equation ((6j). These quantities 
are tabulated in Table [TJ along with other results. For NGC 
2010 (Fig. [T]), for instance, we find!ogt m = 8.12, a m = 0.13, 
and S m = 5%, confirming a focused single SSP solution. 
On the other hand, for NGC 2210 (Fig. 0, Iogt m = 9.74, 
logT m = 10.10, ovrc = 0.79, and S m = 27%, signalling the 
spread in ages discussed above. Section [7] elaborates on how 
these numbers can be used to identify suspicious fits. 
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Figure 1. Left: Observed (black) and best-fit (red) spectra of NGC 2010 for Z = Zq VOOs models. The residual spectrum is shown at 
the bottom. All spectra are in units of the flux at Ao = 4020 A. The "absorption-line" (marked in cyan) at ~ 4205 A is actually a defect 
in the original spectrum, masked from the fits. Top right: Results of the multi-SSP fit. The "bar code" shows the 46 ages in the base. 
Middle right: Best fit Ay and x for single-SSP fits as a function of age. Solid (red) and dashed (orange) vertical lines (which coincide in 
this example) mark the best single-SSP age and mean log i in the multi-SSP fit. Bottom right: rms and A figures of merit as a function 
of age. The best fit model for this Z has an age logt = 8.20 and Ay = 0.00, producing an rms of 0.020 and A = 2.0 percent. Fits with 
log Z/Zq = +0.2 (I.5Z0) yield a miniscule 1% improvement in x 2 ; and logt = 8.11, Ay = 0.11. 
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Figure 2. As Fig.[T] but for NGC 2210 and logZ/Z = -1.7 (0.02Z Q ). Notice that, unlike for NGC 2010, the multi-SSP fit splits into 
2 widely separated ages. The spectral difference between the multi and single population fits is shown as a magenta line in the left panel, 
off-set by -0.1 for clarity. 
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Figure 3. Results of spectral fits of NGC 20f0 with 6 different metallicities. Left: Run of best-fit Ay (bottom), rms (middle) and <5 X 2 
(top) with age. Different lines correspond to different Z's: magenta, blue, green, orange, brown and red correspond to logZ/Zp, = —1.7, 
-1.3, -0.7, -0.4, 0, and +0.2, respectively. Circles mark the best-fit model for each Z. Right: Best fit SSP spectra for all 6 different Z's. 
The observed model is in black and the model and residual spectra are plotted according to the color-scheme above. Note: The Ox — M\ 
residual spectra are multiplied by 3 for clarity. Dashed horizontal lines at y = —0.1, and +0.1 are drawn to facilitate comparisons. The 
best fit model (log Z/Zq = +0.2 in this case) is marked by a thicker window frame. 
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Figure 4. As Fig. |3] but for NGC 2210. 
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5.3 Fits with different metallicities 

The example fits in Figs. [1] and [2] have fixed Z. In Figs.[3]and 
|4]we show how things change for different Z's. The bottom- 
left panels show the best-fit Ay as a function of age for all 6 
metallicities. The middle panel show the rms of the singie- 
SSP spectral fits, while the top panel shows 5 X 2 . Curves are 
color-coded according to Z, and a circle marks the location 
of the best model for the corresponding Z. Panels in the 
right show the best spectral fit achieved for each of the 6 
metallicities. Residual spectra shown in these panels were 
multiplied by 3 for clarity. 

Fig.[3]shows results for NGC 2010. Among all 6 x 46 = 
276 SSPs, the one with log Z/Z Q = +0.2 and logi = 8.11 is 
the one which produces the best fit. Accordingly, this model 
has 8 X 2 = 0. As already noted, the next best Z is solar, for 
which S x 2 — 0.01, i.e., its \ 2 ls j us t 1% worse than the best 
one! Visual inspection of the corresponding spectra confirms 
that these two solutions are indeed indistinguishable. Other 
metallicities produce visibly worse fits, suggesting that Z is 
relatively well constrained to Zq for this SC. 

We note in passing that neither the rms (middle left 
panel) nor A (not plotted) provide a strong numerical dis- 
crimination of different models. As can be seen in Fig. [31 
the best models for different Z's all yield numerically simi- 
lar values for these figures of merit, even when the fits are 
visibly worse. Even if formally equivalent, the S x 2 index is 
much more useful in this sense. 

Results for NGC 2210 are shown in Fig. [3] In this case, 
there is a well defined metallicity: log Z/Zq = —1.7. The 
next best model (log Z/Zq = —1.3) has a nearly twice as 
bad x 2 (^x 2 = 0-94). Fits for higher Z's are so much worse 
that they fall off the top right panel. 



5.4 Age-metallicity-extinction degeneracies 

The plots above already show known couplings between 
parameters. The best way to visualize such covariances is 
through bidimensional A\ 2 maps. 

Fig. [5] shows maps of S x 2 in the three possible projec- 
tions of the t-Z-Av space for NGC 2010. Notice that the 
Z scale is not in physical units. Instead, we prefer to repre- 
sent each of the 6 Z's by a iz = 1 to 6 metallicity index, a 
cosmetic choice which serves to emphasize the coarseness of 
the Z-grid. Dotted contours are for 8 X 2 ^ 1, the inner one 
(S x 2 — 1) being drawn in magenta. The gray-shaded area 
maps the S x 2 ^ 0.5 zone, while the inner 3 contours mark 
the formal 3, 2 and 1 sigma ranges for our choice of Ndof 
(see equation [TJl . 

The well known trade-off between age and metallicity 
is evident in Fig. [5^,. Notice also that solar and over-solar 
models are equally good, as noted before in Fig. [3] where 
best fit spectra for each Z were shown. These two Z's fall 
within the 1 sigma confidence region. In the case of NGC 
2210, the "age-metallicity degeneracy" is not so evident be- 
cause the best-fit metallicity is at the border of the Z-grid 
(Fig. Ho). 

As t and Z assume only the discrete values dictated by 
the set of models used, such Ax 2 (t, Z) maps can be con- 
structed directly from the output of the starlight fits. 
However, starlight only computes the best Ay (and x) for 
each (t, Z) pair, so some extra coding is needed to produce 
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Figure 5. & x 2 contours in the age-metallicity-extinction space 
for NGC 2010. The magenta line correspond to & x 2 = 1 and the 
gray shaded area marks the & x 2 $C 0.5 region. Blue, green and red 
correspond to 3, 2 and 1 sigma confidence regions, respectively. 
Best-fit values are marked with dotted lines. In the top panels a 
dashed line connects the best models for different Z's. 



Testing spectral models for stellar populations 9 




9.5 10 

log t [yr] 




0.5 

A v [mag] 



O 




9.5 - - 



A v [mag] 

Figure 6. As Fig. [5] but for NGC 2210. 



maps involving Ay. We have thus written a complementary 
code which, for each t and Z, forces fits using Ay values 
in a fine grid from to 1. These fits have a single free pa- 
rameter, x, whose optimal value is computed analytically 
from dx 2 /dx = 0. The resulting tables of \ 2 as a function 
of t, Z and Av can then be projected to Ax 2 (t, Av) and 
Ax 2 {Z, A v ) maps. 

As colors redden with increasing age and metallicity, 
one expects negative covariances between Av and both t 
and Z. Such an "age-extinction degeneracy" is indeed ob- 
served in Fig. [5}j (NGC 2010). The effect is also present, but 
much weaker, in NGC 2210 (Fig. IB]:) where the trend is only 
noticeable in the S x 2 0.5 contours. 

For many clusters, contours in the Ay-Z plane show the 
expected "metallicity-extinction degeneracy", but in these 
two examples it is not so clear. A hint of the expected effect is 
present in NGC 2210 (Fig.[6p), visible mainly from the outer 
contours, since the best fit values of Z and Av are at a corner 
of the grid. In NGC 2010, however, the trend is opposite to 
the expected one, with contours showing positive covariance 
(Fig. [SJa). In this case, the blueing of the continuum due 
to the decrease in t as Z increases (Fig. [5^0 is larger than 
the reddening caused by the change in Z, explaining why 
Av increases with Z. This example serves as a reminder 
that the quantitative effects of the combination of t-Z-Ay 
"degeneracies" can be more subtle than those expected on 
simple qualitative grounds. 

These maps show that the Av ^ 1 prior has very little 
effect. This happens because the extinction in our SCs is 
genuinely small. For the same reason, the physical Av 
limit has more effect in constraining Av, as seen, for in- 
stance, in the contours in Fig. (5)]. 

Before closing, we note that, as said in Section f4. 31 the 
statistical interpretation of these surfaces can be adapted to 
other choices for Ndof- For instance, the 5 X 2 — 0.5 shaded 
areas in Figs. [S] and HI] would correspond to 3 sigma contours 
(A^ 2 = 11.8 for a 2D projection) for Ndof — 24, equivalent 
to considering every 34 A of a spectrum as an independent 
datum. 



5.5 Results for all clusters 

Columns 2-5 of Table [1] lists the best t, Z and Av obtained 
for all 27 SCs in our sample. The corresponding value of A is 
listed in column 6. Fig. [7] shows the corresponding spectral 
fits, sorted by the best-fit age. The first cluster shown, NGC 
1818, is actually younger than 100 Myr, the smallest age 
available in the VOOs models, and the fact that the best-fit is 
found for the lowest available Z is an artificial consequence 
of this limitation. Our results for this SC should thus be 
ignored. The next three SCs (NGC 2136, NGC 2214, and 
NGC 1866) are also close to this limit, and results for them 
should be regarded with care, even though the fits are visibly 
better than in the case of NGC 1818. 



6 BAYESIAN ESTIMATES OF STAR 
CLUSTER PROPERTIES 

The fits and A^; 2 maps presented above give an idea of the 
kind of precision in SC parameters achievable with spectral 
fits. So far, however, our assessment of the uncertainties was 
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Figure 7. All spectral fits, sorted by best-fit age. Observed and best fit single SSP spectra are shown in black and red, respectively, with 
the corresponding residual shown in green. Labels to the right of each spectrum list the best values of logt and log Z/Zq. All spectra 
are normalized at Aq = 4020 A, and shifted vertically by f (= one tick mark) from one another. 



mostly qualitative. This section presents a simple Bayesian 
formalism adapted to the problem of SC parameter estima- 
tion. Revised estimates of t, Z and Ay are presented, along 
with their respective uncertainties. 



6.1 Formalism 

For a data set D and a model with parameters p, the pos- 
terior probability of p is given by Bayes theorem: P(p\D) = 
P(p)P(D\p)/P(D). For a P(p) = constant prior and gaus- 
sian errors, P(p\D) oc exp[— ^x 2 (p)L from which estimates 
of any of the individual parameters in p can be obtained 
through marginalization. 



Lets momentarily denote t = logt, Z = log Z/Zq, and 
A = Ay - Hence, p corresponds to (t, Z, A, x), while the data 
D correspond to an observed spectrum (Ox) and its errors 
(e^). The posterior probability distribution function (PDF) 
of, say, the (log) age is then given by 



P(t\D) = J J J P(t,Z,A,x\D)dZdA 



tlx 



(8) 



and similarly for Z, A, and x. 

Computing such multi-dimensional integrals is usually a 
cumbersome task, but with 4 dimensions the problem is still 
manageable and a simple discretization approach suffices. To 
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Table 1. Results for all clusters. Columns 2—5 correspond to the best single SSP fits; columns 6-10 list results of the multi-SSP fits 
(with Z fixed and equal to that given in column 3); columns 11-13 list the Bayesian estimates of t, Z and Av- 



evaluate the PDF of the parameters, one thus needs grids in 
t, Z, A and x. Grids in t and Z are naturally available from 
the SSP models, while grids in A were already computed 
to produce the A\ 2 {t, Av) and A\ 2 (Z, Av) maps shown 
above. The missing a;-grid is the easiest to compute, as the 
effects upon \ 2 °f a change in x with respect to the optimal 
value can be computed analytically. In this computation we 
restrict x to the 0.8-1.2 range, a very generous prior, which, 
however, was not formally imposed in the single-population 
starlight fits presented in Section [5] 

The AZj steps were computed as [Zj+i — Zj-i]/2 for 
intermediate Zj's, while the edge bins (j = 1 and Nz) were 
taken to be symmetric. The same scheme was followed for 
the other parameters, for which the sampling is much finer. 

6.2 Results 

Results for NGC 2010 are shown in Fig. [8] In all panels, the 
solid curve corresponds to PDFs computed with our fiducial 
choice for Ndof (see Section|33J • The best fit value is marked 
by a star. The average value is plotted as a circle, and ±1 
and ±2 sigma error bars (all computed directly from the 
PDF) are also drawn. The correspondance with the contour 
plots in Fig. [5] is evident. The two acceptable metallicities 
translate into the two peaks in P(t). Two peaks are also 
seen in the PDF of Av, but blended. Multi- modal PDFs and 
contour plots like these occur frequently, and are basically a 
by-product of the coarseness of the Z-grid (see Appendix 1X1 



for results obtained with Z- interpolated grids). There are, 
however, exceptions. NGC 2210 is an example whose PDFs 
(Fig. IA2|) are unimodal and centered on the best fit values, 
as can be guessed from Fig. [6] 

Just for illustrative purposes we also draw the PDF 
which would result from adopting a 5 times smaller Ndof ■ As 
expected, the PDFs become broader, but at least in this ex- 
ample, this extreme choice, which tantamounts to equalling 
a full spectral fit to fitting just N\ = 16 observables, the SC 
parameters are still well constrained. 

The standard way of summarizing PDFs is to compute 
their average and standard deviation, as already done in Fig. 
[5] Results for all clusters are listed in the last 3 columns of 
Table[T] These are our official estimates for the SC properties 
and their uncertainties. Note, however, that the table con- 
tains several entries for which <r(log Z) is very small or even 
zero! This is obviously wrong physically and statistically, but 
is a direct mathematical consequence of the coarseness of the 
Z-grid (see Appendix |A")) . In such cases, NGC 2210 being an 
example, the likelihood is so concentrated around one Z that 
neighboring grid points contain essentially no "probability 
mass". When this happens, it is prudent to adopt something 
like half the logZ'-grid separation as a measure of o(\ogZ). 
(The mean half grid separation is 0.2 dex in metallicity and 
0.025 dex in age.) Also, as already noted, results for clus- 
ters too close to the age and metallicity limits of the models 
should be read with caution. 

Typical (sample mean) formal statistical uncertainties 
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Figure 8. Probability distribution functions for t, Z and Av for 
NGC 2010. Solid lines correspond to PDFs with our fiducial choice 
for N^of (i.e., 61), whereas dotted (red) lines are for a 5 times 
smaller Nd j. A star marks the best fit value. One and two sigma 
error bars are also shown, centered on the mean values. In the 
top and middle panels, vertical lines mark ages and metallicities 
available in the set of evolutionary synthesis models used. 



in logf, log Z/Zq and Av are 0.06 dex, 0.09 dex and 0.08 
mag. No trend of uncertainty with mean SC properties was 
found. Notice that the typical uncertainty in logt is nearly 
identical to the 0.05 dex age spacing of the VOOs models, so, 
besides a finer Z-grid, a finer sampling in age would also be 
desirable. 

Simulations which test the ability of the spectral syn- 
thesis methodology used in this paper to recover SC parame- 
ters under different levels of noise are described in Appendix 
[Bl The overall conclusion of these simulations is that the 



method is robust. For S/N > 30 (as most of the data used 
in this paper) one expects bayesian estimates within < 0.1 
dex in age, < 0.2 dex in metallicity and < 0.1 dex in Ay of 
the correct values. 



7 SUMMARY AND DISCUSSION 

We have presented a rigorous but straightforward method- 
ology to fit SC spectra with predictions from modern, pub- 
licly available high resolution evolutionary synthesis mod- 
els. The method is mostly based on the also publicly avail- 
able (but never before used for SCs) spectral synthesis code 
starlight, from which one can derive best t and Z esti- 
mates, as well as likelihood maps in the (t, Z) plane. Objec- 
tive recipes to deal with the absence of errors and account 
for spectral resolution in the definition of the effective num- 
ber of degrees of freedom were presented. A simple bayesian 
approach was followed to compute formal probability distri- 
bution functions for the parameters. The method was ap- 
plied to 27 SCs from the work of Leonardi & Rose (2003), 
modelled with publicly available SSP models of Vazdekis et 
al. (2009) based on the MILES library, "Padova 2000" tracks 
and a Salpeter IMF. 

The main goal of this paper was to setup and illustrate 
this methodology, paving the way for the comparative study 
presented in Paper II, where results obtained with different 
models are presented and compared to independent work. 
This comparison is essential to map systematic uncertain- 
ties, whereas in this paper only statistical uncertainties were 
addressed. 

A first general result is that the spectral fits are very 
good, yielding residuals of the order of A = 2%, equivalent 
to a signal-to-noise ratio of 50. That, of course, does not 
guarantee that the inferred parameters are correct. Insofar 
as Ay is concerned, our results are consistent with previ- 
ous estimates, which indicate very low extinction for SCs in 
this sample (e.g., Bica & Alloin 1986). In Paper II we show 
that, taking CMD determinations of t and Z as a reference, 
spectral fitting does provide reliable results. 

Regarding the estimation of SC properties, the coarse- 
ness of the available models in Z introduces some undesir- 
able complications, like probability islands in the parameter 
space, reflected as multimodal PDFs of age and extinction. 
In cases where one Z is highly preferred, the lack of a finer 
grid leads to an underestimation of the statistical uncertain- 
ties in the SC parameters. Experiments with Z interpolated 
spectra alleviate these problems to some extent, but the ul- 
timate solution is to have a finer sampling in both Z and 
t. 

A curious effect found in some of the multi-SSP fits 
was a split into a light dominating old population (~ 80%) 
and weaker (~ 20%) but significant population of much 
younger age (often the lowest age in the grid), as if the 
SSP model spectra lacked old blue stars. In principle, these 
multi-component starlight fits should not be relevant for 
SCs, but we find that they are able to capture symptoms of 
inadequate single SSP fits in a quantitative manner. Clusters 
with suspicious combinations of different ages may be iden- 
tified by large values of a m (ie., large dispersions in logt). 
Four clusters have a m > 0.6: M15, M79, NGC 2203, and 
NGC 2210 (Figs. [2] and [6j|. M15 and M79 are very similar to 
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Figure 9. Multi-SSP age dispersion plotted in the t-Z plane, 
where (logf) and (log Z/Zq) given by our Baycsian estimates. 
The size of the circles scales with a m (the age dispersion in the 
multi-SSP fits), which ranges from to 0.9 dex. Red circles corre- 
spond to cases where multi-SSP fits are much better than single 
SSP ones (S m > 20%), while green is used when 8 m < 10% and 
blue correspond to cases in between. Large <x m , like those seen 
for leftmost sources (M15, M79 and NGC 2010), is indicative of 
a component not adequately represented in the models, like blue 
HB stars. NGC 1818 is not plotted because it is younger than the 
lower age in the VOOs models. 

NGC 2210. In fact, they all have identical best fit t and Z, 
as well as the same log T m , similar f m , and 5 m values which 
indicate that the multi-SSP fit is substantially better than 
a single SSP fit (5 m = 44% for M15, 97% for M79, and 27% 
for NGC 2210). NGC 2203, on the other hand, is only simi- 
lar to NGC 2210 in terms of a m , and its multi-SSP fit is not 
substantially better than single SSP one (Sm = 8%). Thus, 
the "old plus young" anomaly is only convincingly detected 
in M15, M79 and NGC 2210. 

Following K08, we suggest this result is due to incom- 
plete modelling of the horizontal branch (HB). CMD work 
on M15 (Battistini et al. 1985; Buonanno et al. 1983), M79 
(Kravtsov et al. 1997) and NGC2210 (Brocato et al. 1996) 
have shown that all these clusters have blue HBs, which sup- 
ports our interpretation of the spectral synthesis results. Al- 
ternatively, blue stragglers would have a similar effect. Com- 
binations of UV and optical photometry detect the presence 
of these stars in M15 (Guhathakurta et al. 1993; Ferraro 
& Paresce 1993) and M79 (Lanzoni et al. 2007). The situ- 
ation is less clear in NGC 2210. Based on the similarity of 
NGC 2210 HST CMD to NGC 2257, Johnson et al. (1999) 
suggest that the cluster have both a blue HB and blue strag- 
glers, while Brocato et al (1996) have argued (on the basis 
of a ground based CMD) that the stars in the blue straggler 
region are likely field contaminants in NGC 2210. 

In broad terms, one expects bluer HB at lower Z (Harris 
1996; Greggio & Renzini 1990), and hence the age-spread 
noticed in the multi-SSP fits should increase for decreasing 
Z. Fig. [9] presents tentative evidence of this relation. SCs 
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Figure 10. Extension of statistically indistinguishable fits to the 
3650-4600 A data to the red limit of the models for 4 SCs. For 
each cluster, the red and green lines are fits with different metal- 
licities (log Z/Z Q = -1.3 and -0.7 for NGC 1651 and NGC 1795, 
and log Z/Z Q = and +0.2 for NGC 2133 and NGC 2010). The 
observed spectra are plotted in black. 

with large <j m and S m tend to be those with smallest Z. 
Because Z and t are strongly anti-correlated in this sample, 
the highest a m clusters are also among the oldest ones. Given 
that Z is not the single factor determining HB morphology, 
not to mention other possibilities (like blue stragglers and 
stochastic effects) and uncertainties involved, we shall not 
elaborate further on this issue. 

Moving away from caveats in the stellar population is- 
sues, we must acknowledge at least one limitation of our 
results: Wavelength coverage. The spectra used here cover 
a relatively small region in the blue-near-UV. Not only it 
misses some strong t and Z sensitive absorption features, 
like H/3 and the Mgl lines, but the lack of a wider baseline 
prevents a firmer handle on Ay - No doubt a wider coverage 
would help better constrain SC properties, but how much? 
Could more data help, for instance, lifting the ambiguities 
in Z found for several SCs? 

A qualitative answer to this question is given in Fig. 1101 
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where, for each of four SCs, we plot the two best SSP fits 
with different Z, extending it from the 3650-4600 A range 
covered by the data to the red limit of the VOOs models. For 
the two clusters at the bottom of the plot, NGC 1651 and 
NGC 1795, for which log Z/Z Q = -1.3 and -0.7 provide 
statistically indistinguishable fits, the extrapolated model 
spectra differ a lot both in continuum shape and in the depth 
of absorption features to the red of 4600 A. In these cases, 
data in the red would clearly help distinguishing the fits. 
For NGC 2010 and NGC 2133, on the other hand, spectra 
in the red would not be as useful, as the extrapolated mod- 
els are essentially indistinguishable. A wider baseline would 
be needed to pin down the best solution in these younger 
systems. 
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Figure Al. Mean percentage spectral residual (A) as a function 
of age for Z-interpolated models as compared to the original ones 
for log Z/Z Q = —1.3, -0.7, -0.4 and 0.0. Average values of A for 
all ages arc listed in the inset table. The message from this plot is 
that interpolating SSP spectra in Z introduces errors which are of 
the same order as the spectral residuals obtained in fits of actual 
SC data. 



APPENDIX A: INTERPOLATIONS IN 
METALLICITY 

In the absence of a fine Z-grid, one may resort to inter- 
polations. Because of stellar evolution effects, interpolating 
in Z at fixed t is not strictly valid, but one expects it to 
work at least approximately. This appendix describes re- 
sults of experiments with Z-interpolated models. The orig- 
inal Nz — 6 grid was interpolated to Nz = 39 with 
log Z/Z between -1.7 and +0.2 in 0.05 dex steps. The 
spectral interpolations were carried out in log-log space, with 
log F x {t;Z) = ax(t)logZ/Z Q +bx(t). 

To test how well Z-interpolation works we have com- 
pared the original log Z/Z Q = —1.3, —0.7, —0.4 and 0.0 
spectra with those obtained through interpolation of mod- 
els with adjacent Z- values (for instance, using the —1.3 and 
—0.4 models to interpolate to —0.7). The results are shown 
in Fig. IA1I where we plot the mean percentage spectral de- 
viation A against t for each of these four intermediate Z's. 
The residuals are computed after re-normalizing the spectra 
at 4020 A, eliminating differences in the absolute flux scale. 
The plot shows that interpolation leads to typical errors of 
1% in terms of spectral residuals, but that in many cases A 
approaches or even exceeds 2%. These residuals are of the 
same order of the A values obtained in fits of actual SCs 
(Table [TJ ! Naturally, in these tests the interpolations are 
based on models farther apart in Z than when using the full 
grid. Nevertheless, we take these results as a signal that an 
interpolated grid may not lead to as much improvement as 
one would hope to achieve. 

We now check how a finer Z-grid changes our parame- 
ter estimates. Fig. IA2l shows the marginalized PDFs of t, Z 
and A v for NGC 2010, NGC 2210 and NGC 1651. The plot 
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Figure A2. Comparison of the t, Z and Ay PDFs obtained 
with the original (dotted red lines) and Z-interpolated grids (solid 
black), for three SCs in our sample. As in Fig. \E\ a star marks 
the best-fit value, and the filled circle marks bayesian estimate. 
Error bars are also shown. 



compares the PDFs derived with the original Nz = 6 Z-grid 
(dotted red lines) with those derived with the Nz — 39 inter- 
polated models (solid black). For NGC 2010 and NGC 2210 
the Z-interpolated PDFs do not differ significantly from the 
original ones. The clearest changes are in the PDF(Av) for 
NGC 2010, which becomes smoother (single-peaked), and 
in the PDF(Z) for NGC 2210, which rises steeply towards 
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log Z/Zq — — 1.7 (confirming that the lowest Z models are 
indeed the ones which better fit the data for this SC). More 
significant changes occur in the case of NGC 1651. With the 
Z-interpolated models, PDF(£) has a pronounced peak in 
between the two peaks obtained with the original grid. Sim- 
ilarly, the PDF(i?) now peaks in between the two original Z 
values with highest probability. The more dramatic change 
is in PDF(Av), whose original two peaks merge into a single 
one at an intermediate value of Ay. For all three parameters 
the Z-interpolated models produce more focused PDFs, and 
thus smaller formal uncertainties (see error bars in the NGC 
1651 panels of Fig. IA2[l , Notwithstanding these apparently 
large differences, the bayesian estimates of t, Z and Ay de- 
rived with the coarser original Z-grid are consistent with 
those derived with the better sampled interpolated Z-grid. 
They are also more conservative, in the sense of having larger 
error bars. 

Considering the sample as a whole, we detect no bias 
between parameters estimated with the original and Z- 
interpolated grids. The rms difference in the bayesian es- 
timates of log£ and logZ/^© are just 0.08 and 0.12 dex 
respectively, while for Ay the rms difference is 0.08 mag. 
Z-interpolated grids produce better spectral fits, but only 
marginally so. In terms of the A figure of merit, the Z- 
interpolated fits yield a sample-average value of A = 2.39%, 
compared to 2.44% obtained with the original grid. 

Regarding formal (PDF-based) parameter uncertain- 
ties, our experiments show that when an uncertainty de- 
rived with the original grid is "large" (> 0.05 dex in log£, 

> 0.15 dex in logZ/Z©, and > 0.1 mag in Ay), it tends 
to be smaller with the Z-interpolated grid, and vice versa. 
This behavior is consistent with qualitative expectations. In 
cases where a PDF is unresolved due to a coarse grid, lead- 
ing to a formal error of zero, interpolation resolves the PDF, 
producing a measurable second moment (standard deviation 

> 0). For instance, for NGC 2210 we find cr(logZ) = with 
the original grid (all probability is in the first Z-bin), and 
0.02 dex with the interpolated one. Conversely, when, say, 
the PDF of Z is spread over two or more adjacent origi- 
nal Z-bins, a finer grid may produce a more focused (nar- 
rower) PDF, as in the case of NGC 1651. In any case, these 
are relatively small differences, such that the sample-average 
uncertainties derived with both methods are virtually indis- 
tinguishable. 

All in all, we see no major gain in resorting to Z- 
interpolated models. Further considering the uncertainties 
introduced by such interpolations (Fig. IA1[) . we think it is 
safer to wait for properly computed, finer Z-grids to be re- 
leased by than to interpolate in Z. 



APPENDIX B: EFFECTS OF DATA QUALITY 

The starlight code, used throughout this paper, has been 
extensively used and tested in the context of mixed stellar 
populations (galaxies), but this is the first time it is used 
to fit SCs. This appendix presents simulations designed to 
evaluate the reliability of SC properties derived with this 
code under different levels of noise. 

To test the effects of data quality, we have performed 
simulations in which white gaussian noise is added to a 
model SSP spectrum of known t, Z and Ay. These per- 



turbed spectra are processed exactly as the observed ones, 
and input and output parameters are then compared. The 
input models were built using the VOOs SSPs which best 
fit the 27 actual SCs in our sample, whose properties are 
listed in columns 2-4 of Table [1] Their corresponding spec- 
tra were taken directly from starlight's outpu10, and thus 
mimic the actual data in terms of wavelength coverage and 
resolution. Five realizations of a noise spectrum (n\) were 
generated, and five others were derived from these by in- 
verting the sign of n\. These 10 perturbation spectra were 
scaled to emulate S/N ratios between 10 and 95. 

Fig. IB1I illustrates the results for three test SCs, with 
parameters based on those found for NGC 1651, NGC 1795 
and NGC 2010. All panels show output minus input param- 
eters against the S/N. The left and middle columns show, 
for each S/N, the results of an analysis based on an indi- 
vidual spectrum. The examples picked for these panels differ 
only in the sign of the perturbation, such that if in the left 
column the spectrum feed into starlight is F\ + n\, in the 
middle one the analysis is carried out with F\ — n\. Also, 
n\ has the same shape for all S/N values, such that the 
only thing that changes along the x-axis is the amplitude 
of the noise. This cosmetic choice is done only for presen- 
tation purposes, as it produces more illustrative plots than 
obtained picking perturbed spectra at random. Finally, the 
third column of plots shows statistical results considering all 
10 spectra available at each S/N . 

Consider first the left and middle panels. Stars show the 
difference between best fit and input values of log t, log Z 
and Ay as a function of S/N . The solid line shows the dif- 
ferences between the bayesian estimates and the input val- 
ues, while the grey shaded areas show the corresponding ± 
one sigma formal (PDF-based) uncertainty. Besides the ex- 
pected convergence as S/N increases, these plots illustrate 
covariances between age, metallicity and extinction. When 
the ^-estimate does not change, underestimated t's are com- 
pensated by overestimated Ay's, and vice- versa, as in the 
left panels for NGC 1651 and in the left and middle panels 
for NGC 1795. In the middle panel for NGC 1651, one sees 
that both t and Ay are overestimated at low S/N due to 
the fact that fits with the logZ/Z® = —1.3 grid point are 
preferred over the logZ/Z® — —0.7 original value. In this 
particular example, only for S/N > 30 the input parameters 
are well recovered. In the case of NGC 2010, the left panels 
show that solutions oscillate between Z = 1.5Zq (the input 
value) and Zq, and all SC parameters vary in concert. This 
metallicity confusion occurs for S/N as large as 50, which is 
not surprising, given the high degree of spectral similarity 
between solar and over-solar VOOs models in this age range 
(see Fig. [3}. In the example examined in the middle panels, 
Z comes out badly underestimated for S/N ^ 15. Above 
that Z is well recovered, but minor differences in Ay and t 
can persist to S/N as high as 70. 

These examples illustrate what may happen in the anal- 
ysis of individual spectra observed under different regimes of 
S/N . The general message spelt by these simulations is that, 
when dealing with one single object and using only the kind 
of observational data used in this paper, one should be care- 
ful about biases in parameter estimates obtained through 



2 Only version 5 of the code outputs the best fit SSP spectrum. 
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Figure Bl. Illustration of the effects of noise upon derived SC 
properties for simulations based on the t, Z and Ay values ob- 
tained for NGC 1651 (top panels), NGC 1795 (middle) and NGC 
2010 (bottom). The 1st and 2nd columns of plots show results of 
the analysis of two individual perturbed spectra. Error bars and 
the shaded area show the PDF based logt ± o-(logf), and simi- 
larly for log Z and Ay, while the (green) star marks the best-fit 
value. The last column shows the statistics of the simulations. In 
this case the shaded area shows the mean and ± one sigma range 
of the PDF based values of logt, logZ and Ay. The points with 
error bars (in red) show the mean and standard deviation of the 
best-fit values among 10 realizations of the noise for each S/N. 



spectral synthesis if S/N < 30. Beyond this S/N, output 
minus input deviations are too small to worry about. 

The right panels in Fig. IB II show the results of a statis- 
tical analysis of all perturbed spectra for each test-SC. The 
yellow shaded ranges marks the dispersion in the bayesian 
estimates of logt, \ogZ, and Ay among the 10 perturbed 
versions of a same spectrum. Similarly, the red error bars 
show the dispersion in the best fit estimates as a function 
of S/N. The sometimes large deviations observed in the left 
and middle panels average out in these statistical results. 
Biases in the estimates are still present at low S/N, but al- 
most always within one sigma of output = input. The most 
noticeable bias in the right panels of Fig. IB II occurs for Ay 
in NGC 1651, and only for low S/N. Since Ay < was for- 
bidden in the fits, and the input value of Ay is exactly zero 
in this case, one can only err upwards, so this bias can be 
understood as an edge effect. 

Fig. IB2I summarizes the statistics of the noise-induced 
variations in bayesian parameter estimates for all 27 test 
SCs. The top panels, show the standard deviation (consid- 
ering the 10 perturbations) of the log i estimate against the 
input value of log t for three values of S/N: 10 (red open cir- 
cles), 20 (green open squares) and 30 (blue filled triangles). 
Dotted lines mark the average a (logt) for these 3 values of 
S/N. The plot shows that the typical uncertainties in logf 
due to noise are of < 0.2 dex for S/N > 10. For SCs in 
the 9.0 to 9.5 logt range the uncertainty reaches 0.3 dex for 
S/N — 10, which is acceptable considering such low quality 
data. The middle and bottom panels show the results for 
log Z/Zq and Ay. For the metallicity, noise at a S/N = 10 
level introduces up to 0.5 dex uncertainty, but the overall 
average value is 0.3 dex. As for the ages, the uncertainty in 
Z peaks at intermediate values. A plausible explanation for 
this is that in the middle of a (t or Z) grid one can look 
to either side and find reasonably similar spectra, whereas 
closer to the edges there is not so much option, thus reduc- 
ing the ranges in t and Z where statistically acceptable fits 
can be found. For Ay we find up to 0.2 mag uncertainty, 
with a typical value of 0.1 mag. 

The overall conclusion of these experiments is that, for 
the kind of data used in this paper (medium resolution spec- 
tra in the blue range), S/N > 30 is required to accurate SC 
parameter estimates from spectral synthesis of an individual 
object. When large samples of objects are analyzed, one can 
trust statistical results obtained with lower quality data. 

We emphasize that our conclusion that S/N > 30 is re- 
quired to produce accurate parameter estimates is not equiv- 
alent to saying the method fails at lower S/N. Our simula- 
tions show that, as expected, output minus input differences 
in SC properties tend to increase as S/N decreases, but so 
do the associated error bars. The reliability of a parame- 
ter estimation method is better measured not in terms of 
absolute differences, but in terms of differences in units of 
the associated uncertainty, e.g., [logt — rogt(input)]/<j(logt). 
Evaluated in this way, spectral synthesis is capable of pro- 
ducing good results for any of the S/N investigated. This 
is illustrated in Fig. IB3I which shows histograms of output 
minus input parameters in units of the corresponding formal 
uncertainty. All 10 perturbed versions of all 27 SCs are in- 
cluded in these histograms. Bottom, middle and top panels 
show results for S/N = 10, 20 and 30, respectively. Dot- 
ted lines draw gaussians of unit a and normalized to the 
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Figure B3. Distribution of the difference between the output 
bayesian estimate and the input parameter in units of the formal 
uncertainty for logt (left), logZ middle, and Ay (right). Vertical 
dashed lines mark the average value of the normalized difference. 
Dotted lines (in red) correspond to gaussians of unit standard 
deviation and normalized to the same area as the simulated data. 
The histograms include results for 270 runs, corresponding to 10 
perturbations of the 27 test-SCs. Botttom, middle and top panels 
correspond to simulations with S/N = 10, 20 and 30, respectively. 



pected. This happens because as the noise decreases the 
PDF of Z becomes increasingly focused and centered on the 
correct grid- value of Z, and hence \ogZ = log Z(input). 



Figure B2. Statistical uncertainties on log t (top panel), log.Z 
(middle) and Ay (bottom) for simulations with S/N = 10 (red 
open circles), 20 (green open squares) and 30 (blue solid trian- 
gles). Dotted horizontal lines mark the mean value over all 27 sim- 
ulated clusters. Each point corresponds to results from 10 Monte 
Carlo realizations of the noise. The 27 points (for each S/N level) 
correspond to test clusters whose parameters are identical to those 
in columns 2-4 of Table [T] 



same area. The histograms follow quite well the the gaus- 
sian expected in the ideal case, with the noticeable exception 
of Ay, which is systematically overestimated as a result of 
the already mentioned Ay edge effect. For S/N — 10 
this bias is close to 1 sigma, as can be seen by the loca- 
tion of the vertical dashed line, which marks the average 
value of {Ay~ - A v (input)] /cr(A v ). For S/N = 20 this dif- 
ference reduced to 0.5 sigma on average, while for better 
data the distribution of [Ay — yly(input)]/cr(Av) becomes 
progressively less skewed. Another effect which can be no- 
ticed in this plot is that for high S/N the distribution of 
[logZ — log Z(input)]/<j(log Z) becomes narrower than ex- 



